#' ---
#' title: "Main Analysis with Ordered Logit"
#' author: "Gento Kato & Fan Lu"
#' date: "June 8, 2023"
#' ---
#' 
#' 
#' # Preparation 
#' 

## Clean Up Space
rm(list=ls())

## Set Working Directory (Automatically) ##
setwd(dirname(rstudioapi::getActiveDocumentContext()$path)); 

require(interflex)
require(estimatr)
require(texreg)
require(ggplot2)
require(ggrepel)
require(haven)
require(MASS)
require(sandwich)
require(lmtest)

survey22 <- readRDS("data_survey22_v7.rds")
icpp09 <- readRDS("data_icpp09_v7.rds")

survey22$immigincrease <- factor(survey22$immigincrease, 
                               labels = c("Disagree","Somewhat Disagree",
                                          "Somewhat Agree", "Agree"))
survey22$foreignsuff <- factor(survey22$foreignsuff, 
                               labels = c("Don't Think So",
                                          "Don't Think Somewhat",
                                          "Can't Say",
                                          "Think Somewhat",
                                          "Think So"))
icpp09$immigincrease <- factor(icpp09$immigincrease, 
                               labels = c("Disagree","Somewhat Disagree",
                                          "Somewhat Agree", "Agree"))
icpp09$foreignsuff <- factor(icpp09$foreignsuff, 
                             labels = c("Don't Think So",
                                        "Don't Think Somewhat",
                                        "Can't Say",
                                        "Think Somewhat",
                                        "Think So"))

cfmap <- list("(Intercept)"="(Intercept)",
              "edu2"="University Education",
              "cohortCohort II (18+ in 1976-1989)"="Cohort II (18+ in 1975-89)",
              "cohortCohort III (18+ in 1990-99)"="Cohort III (18+ in 1990-99)",
              "cohortCohort IV (18+ in 2000-09)"="Cohort IV (18+ in 2000-09)",
              "cohortCohort V (18+ in 2010-)"="Cohort V (18+ in 2010-)",
              "edu2:cohortCohort II (18+ in 1976-1989)"="University * Cohort II",
              "edu2:cohortCohort III (18+ in 1990-99)"="University * Cohort III",
              "edu2:cohortCohort IV (18+ in 2000-09)"="University * Cohort IV",
              "edu2:cohortCohort V (18+ in 2010-)"="University * Cohort V",
              "male"="Gender (Male)", 
              "cohortCohort II (18+ in 1976-1989):male"="Male * Cohort II",
              "cohortCohort III (18+ in 1990-99):male"="Male * Cohort III",
              "cohortCohort IV (18+ in 2000-09):male"="Male * Cohort IV",
              "cohortCohort V (18+ in 2010-):male"="Male * Cohort V",
              "incomecatMiddle" = "Income (Middle)",
              "cohortCohort II (18+ in 1976-1989):incomecatMiddle"="Income (Middle) * Cohort II",
              "cohortCohort III (18+ in 1990-99):incomecatMiddle"="Income (Middle) * Cohort III",
              "cohortCohort IV (18+ in 2000-09):incomecatMiddle"="Income (Middle) * Cohort IV",
              "cohortCohort V (18+ in 2010-):incomecatMiddle"="Income (Middle) * Cohort V",
              "incomecatHigh" = "Income (High)",
              "cohortCohort II (18+ in 1976-1989):incomecatHigh"="Income (High) * Cohort II",
              "cohortCohort III (18+ in 1990-99):incomecatHigh"="Income (High) * Cohort III",
              "cohortCohort IV (18+ in 2000-09):incomecatHigh"="Income (High) * Cohort IV",
              "cohortCohort V (18+ in 2010-):incomecatHigh"="Income (High) * Cohort V",
              "incomecatMissing" = "Income (Missing)",
              "cohortCohort II (18+ in 1976-1989):incomecatMissing"="Income (Missing) * Cohort II",
              "cohortCohort III (18+ in 1990-99):incomecatMissing"="Income (Missing) * Cohort III",
              "cohortCohort IV (18+ in 2000-09):incomecatMissing"="Income (Missing) * Cohort IV",
              "cohortCohort V (18+ in 2010-):incomecatMissing"="Income (Missing) * Cohort V",
              "workstatStudent/Housemaker/Part-Time" = 
                "Student/Housemaker/Part-Time",
              "cohortCohort II (18+ in 1976-1989):workstatStudent/Housemaker/Part-Time"=
                "Student/Housemaker/Part-Time * Cohort II",
              "cohortCohort III (18+ in 1990-99):workstatStudent/Housemaker/Part-Time"=
                "Student/Housemaker/Part-Time * Cohort III",
              "cohortCohort IV (18+ in 2000-09):workstatStudent/Housemaker/Part-Time"=
                "Student/Housemaker/Part-Time * Cohort IV",
              "cohortCohort V (18+ in 2010-):workstatStudent/Housemaker/Part-Time"=
                "Student/Housemaker/Part-Time * Cohort V",
              "workstatStudent/Part-Time" = 
                "Student/Part-Time",
              "cohortCohort II (18+ in 1976-1989):workstatStudent/Part-Time"=
                "Student/Part-Time * Cohort II",
              "cohortCohort III (18+ in 1990-99):workstatStudent/Part-Time"=
                "Student/Part-Time * Cohort III",
              "cohortCohort IV (18+ in 2000-09):workstatStudent/Part-Time"=
                "Student/Part-Time * Cohort IV",
              "cohortCohort V (18+ in 2010-):workstatStudent/Part-Time"=
                "Student/Part-Time * Cohort V",
              "workstatSelf-Employed/Full-Time/Managerial" = 
                "Self-Employed/Full-Time",
              "cohortCohort II (18+ in 1976-1989):workstatSelf-Employed/Full-Time/Managerial"=
                "Self-Employed/Full-Time * Cohort II",
              "cohortCohort III (18+ in 1990-99):workstatSelf-Employed/Full-Time/Managerial"=
                "Self-Employed/Full-Time * Cohort III",
              "cohortCohort IV (18+ in 2000-09):workstatSelf-Employed/Full-Time/Managerial"=
                "Self-Employed/Full-Time * Cohort IV",
              "cohortCohort V (18+ in 2010-):workstatSelf-Employed/Full-Time/Managerial"=
                "Self-Employed/Full-Time * Cohort V",
              "married"="Currently Married",
              "cohortCohort II (18+ in 1976-1989):married"="Married * Cohort II",
              "cohortCohort III (18+ in 1990-99):married"="Married * Cohort III",
              "cohortCohort IV (18+ in 2000-09):married"="Married * Cohort IV",
              "cohortCohort V (18+ in 2010-):married"="Married * Cohort V",
              "urban"="Urbanness (Current Residence)",
              "cohortCohort II (18+ in 1976-1989):urban"="Urbanness * Cohort II",
              "cohortCohort III (18+ in 1990-99):urban"="Urbanness * Cohort III",
              "cohortCohort IV (18+ in 2000-09):urban"="Urbanness * Cohort IV",
              "cohortCohort V (18+ in 2010-):urban"="Urbanness * Cohort V")

#'
#' # Analysis
#'

#+
#########################################
## Analysis of 2009 Survey (Figure A5) ##
#########################################

## Models with Categorical Interaction
mA_icpp09 <- polr(immigincrease~edu2*cohort+
                         male*cohort+
                         incomecat*cohort+
                         workstat*cohort+
                         married*cohort+
                         urban*cohort, 
                       data=icpp09)
summary(mA_icpp09)
out_mA_icpp09 <- coeftest(mA_icpp09, 
                          vcov. = vcovCL(mA_icpp09, cluster = ~ area))
ci_mA_icpp09 <- coefci(mA_icpp09, 
                          vcov. = vcovCL(mA_icpp09, cluster = ~ area))
mB_icpp09 <- polr(foreignsuff~edu2*cohort+
                    male*cohort+
                    incomecat*cohort+
                    workstat*cohort+
                    married*cohort+
                    urban*cohort, 
                  data=icpp09)
out_mB_icpp09 <- coeftest(mB_icpp09, 
                          vcov. = vcovCL(mB_icpp09, cluster = ~ area))
ci_mB_icpp09 <- coefci(mB_icpp09, 
                       vcov. = vcovCL(mB_icpp09, cluster = ~ area))

screenreg(list(mA_icpp09,mB_icpp09),
          override.se = list(out_mA_icpp09[,2],out_mB_icpp09[,2]),
          override.pvalues = list(out_mA_icpp09[,4],out_mB_icpp09[,4]),
          digits=3, include.ci=FALSE, single.row=TRUE,
          custom.model.names = c("Admit More Foreigners","Integrate Foreign Residents"),
          custom.coef.map = cfmap)

## For Plotting ##

mA_icpp09_list <- list()
out_mA_icpp09_list <- list()
ci_mA_icpp09_list <- list()
mB_icpp09_list <- list()
out_mB_icpp09_list <- list()
ci_mB_icpp09_list <- list()

for(i in 1:4) {

  mA_icpp09_list[[i]] <- polr(immigincrease~edu2*relevel(cohort,levels(cohort)[i])+
                      male*relevel(cohort,levels(cohort)[i])+
                      incomecat*relevel(cohort,levels(cohort)[i])+
                      workstat*relevel(cohort,levels(cohort)[i])+
                      married*relevel(cohort,levels(cohort)[i])+
                      urban*relevel(cohort,levels(cohort)[i]), 
                    data=icpp09)
  out_mA_icpp09_list[[i]] <- coeftest(mA_icpp09_list[[i]], 
                            vcov. = vcovCL(mA_icpp09_list[[i]], cluster = ~ area))
  ci_mA_icpp09_list[[i]] <- coefci(mA_icpp09_list[[i]], 
                         vcov. = vcovCL(mA_icpp09_list[[i]], cluster = ~ area))
  mB_icpp09_list[[i]] <- polr(foreignsuff~edu2*relevel(cohort,levels(cohort)[i])+
                      male*relevel(cohort,levels(cohort)[i])+
                      incomecat*relevel(cohort,levels(cohort)[i])+
                      workstat*relevel(cohort,levels(cohort)[i])+
                      married*relevel(cohort,levels(cohort)[i])+
                      urban*relevel(cohort,levels(cohort)[i]), 
                    data=icpp09)
  out_mB_icpp09_list[[i]] <- coeftest(mB_icpp09_list[[i]], 
                            vcov. = vcovCL(mB_icpp09_list[[i]], cluster = ~ area))
  ci_mB_icpp09_list[[i]] <- coefci(mB_icpp09_list[[i]], 
                         vcov. = vcovCL(mB_icpp09_list[[i]], cluster = ~ area))

}

mA_edu2_icpp09 <- 
  as.data.frame(rbind(
    c(out_mA_icpp09_list[[1]][1,c(1,4)],ci_mA_icpp09_list[[1]][1,]),
    c(out_mA_icpp09_list[[2]][1,c(1,4)],ci_mA_icpp09_list[[2]][1,]),
    c(out_mA_icpp09_list[[3]][1,c(1,4)],ci_mA_icpp09_list[[3]][1,]),
    c(out_mA_icpp09_list[[4]][1,c(1,4)],ci_mA_icpp09_list[[4]][1,])
  ))
colnames(mA_edu2_icpp09) <- c("est","p","lci","uci")
mA_edu2_icpp09$dv <- "Admit More Foreigners"
# unique(sort(icpp09$univyr))
# median(1948:1974); median(1975:1989); median(1990:1999); median(2000:2008) 
mA_edu2_icpp09$x <- c(1961,1982,1994.5,2004)

mB_edu2_icpp09 <- 
  as.data.frame(rbind(
    c(out_mB_icpp09_list[[1]][1,c(1,4)],ci_mB_icpp09_list[[1]][1,]),
    c(out_mB_icpp09_list[[2]][1,c(1,4)],ci_mB_icpp09_list[[2]][1,]),
    c(out_mB_icpp09_list[[3]][1,c(1,4)],ci_mB_icpp09_list[[3]][1,]),
    c(out_mB_icpp09_list[[4]][1,c(1,4)],ci_mB_icpp09_list[[4]][1,])
  ))
colnames(mB_edu2_icpp09) <- c("est","p","lci","uci")
mB_edu2_icpp09$dv <- "Integrate Foreign Residents"
mB_edu2_icpp09$x <- c(1961,1982,1994.5,2004)

m_edu2_icpp09 <- rbind(mA_edu2_icpp09,mB_edu2_icpp09)
m_edu2_icpp09$dv <- factor(m_edu2_icpp09$dv,levels=unique(m_edu2_icpp09$dv))
m_edu2_icpp09$lb <- factor(c("I","II","III","IV"), levels=c("I","II","III","IV"))
m_edu2_icpp09

#+ fig.width=7, fig.height=4
## Figure
ggplot(m_edu2_icpp09) + 
  geom_vline(aes(xintercept=1975), linetype=1, color="gray80") + 
  geom_vline(aes(xintercept=1990), linetype=1, color="gray80") + 
  geom_vline(aes(xintercept=2000), linetype=1, color="gray80") + 
  geom_vline(aes(xintercept=2010), linetype=1, color="gray80") + 
  geom_hline(aes(yintercept=0), linetype=2, color="gray50") + 
  geom_errorbar(aes(x=x, ymin=lci,ymax=uci), alpha=0.8, width=2) + 
  geom_point(aes(x=x, y=est), size=1.5, color="gray10", alpha=0.8) + 
  geom_text(aes(x=x, y=-0.7, label=lb), family="serif") + 
  facet_grid(.~dv) + 
  scale_fill_manual(values=c("black","gray70")) + 
  scale_x_continuous(breaks=c(1975,1990,2000,2010),limits = c(1948,2010)) +
  # coord_cartesian(ylim=c(-0.02,0.3)) + 
  labs(x = "Year of Turning 19 (Typical Year of Entering University)",
       y = "Conditional Coefficient with 95% Confidence Interval") +
  theme_classic(base_family="serif") + 
  theme(plot.title = element_text(hjust=0.5),
        panel.background = element_rect(color="black", size=1),
        strip.text = element_text(size=11, face="bold"),
        legend.position="none")
ggsave("mx_cplot_icpp09_ologit.pdf", width=7, height=4)

#+
#########################################
## Analysis of 2022 Survey (Figure A6) ##
#########################################

## Models with Categorical Interaction
mA_survey22 <- polr(immigincrease~edu2*cohort+
                    male*cohort+
                    incomecat*cohort+
                    workstat*cohort+
                    married*cohort+
                    urban*cohort, 
                  data=survey22)
out_mA_survey22 <- coeftest(mA_survey22, 
                          vcov. = vcovCL(mA_survey22))
ci_mA_survey22 <- coefci(mA_survey22, 
                       vcov. = vcovCL(mA_survey22))
mB_survey22 <- polr(foreignsuff~edu2*cohort+
                    male*cohort+
                    incomecat*cohort+
                    workstat*cohort+
                    married*cohort+
                    urban*cohort, 
                  data=survey22)
out_mB_survey22 <- coeftest(mB_survey22, 
                          vcov. = vcovCL(mB_survey22))
ci_mB_survey22 <- coefci(mB_survey22, 
                       vcov. = vcovCL(mB_survey22))

screenreg(list(mA_survey22,mB_survey22),
          override.se = list(out_mA_survey22[,2],out_mB_survey22[,2]),
          override.pvalues = list(out_mA_survey22[,4],out_mB_survey22[,4]),
          digits=3, include.ci=FALSE, single.row=TRUE,
          custom.model.names = c("Admit More Foreigners","Integrate Foreign Residents"),
          custom.coef.map = cfmap)

## For Plotting ##

mA_survey22_list <- list()
out_mA_survey22_list <- list()
ci_mA_survey22_list <- list()
mB_survey22_list <- list()
out_mB_survey22_list <- list()
ci_mB_survey22_list <- list()

for(i in 1:4) {
  
  mA_survey22_list[[i]] <- polr(immigincrease~edu2*relevel(cohort,levels(cohort)[i])+
                                male*relevel(cohort,levels(cohort)[i])+
                                incomecat*relevel(cohort,levels(cohort)[i])+
                                workstat*relevel(cohort,levels(cohort)[i])+
                                married*relevel(cohort,levels(cohort)[i])+
                                urban*relevel(cohort,levels(cohort)[i]), 
                              data=survey22)
  out_mA_survey22_list[[i]] <- coeftest(mA_survey22_list[[i]], 
                                      vcov. = vcovCL(mA_survey22_list[[i]]))
  ci_mA_survey22_list[[i]] <- coefci(mA_survey22_list[[i]], 
                                   vcov. = vcovCL(mA_survey22_list[[i]]))
  mB_survey22_list[[i]] <- polr(foreignsuff~edu2*relevel(cohort,levels(cohort)[i])+
                                male*relevel(cohort,levels(cohort)[i])+
                                incomecat*relevel(cohort,levels(cohort)[i])+
                                workstat*relevel(cohort,levels(cohort)[i])+
                                married*relevel(cohort,levels(cohort)[i])+
                                urban*relevel(cohort,levels(cohort)[i]), 
                              data=survey22)
  out_mB_survey22_list[[i]] <- coeftest(mB_survey22_list[[i]], 
                                      vcov. = vcovCL(mB_survey22_list[[i]]))
  ci_mB_survey22_list[[i]] <- coefci(mB_survey22_list[[i]], 
                                   vcov. = vcovCL(mB_survey22_list[[i]]))
  
}

mA_edu2_survey22 <- 
  as.data.frame(rbind(
    c(out_mA_survey22_list[[1]][1,c(1,4)],ci_mA_survey22_list[[1]][1,]),
    c(out_mA_survey22_list[[2]][1,c(1,4)],ci_mA_survey22_list[[2]][1,]),
    c(out_mA_survey22_list[[3]][1,c(1,4)],ci_mA_survey22_list[[3]][1,]),
    c(out_mA_survey22_list[[4]][1,c(1,4)],ci_mA_survey22_list[[4]][1,])
  ))
colnames(mA_edu2_survey22) <- c("est","p","lci","uci")
mA_edu2_survey22$dv <- "Admit More Foreigners"
# unique(sort(survey22$univyr))
# median(1963:1989); median(1990:1999); median(2000:2009); median(2010:2022) 
mA_edu2_survey22$x <- c(1976,1994.5,2004.5,2016)

mB_edu2_survey22 <- 
  as.data.frame(rbind(
    c(out_mB_survey22_list[[1]][1,c(1,4)],ci_mB_survey22_list[[1]][1,]),
    c(out_mB_survey22_list[[2]][1,c(1,4)],ci_mB_survey22_list[[2]][1,]),
    c(out_mB_survey22_list[[3]][1,c(1,4)],ci_mB_survey22_list[[3]][1,]),
    c(out_mB_survey22_list[[4]][1,c(1,4)],ci_mB_survey22_list[[4]][1,])
  ))
colnames(mB_edu2_survey22) <- c("est","p","lci","uci")
mB_edu2_survey22$dv <- "Integrate Foreign Residents"
mB_edu2_survey22$x <- c(1976,1994.5,2004.5,2016)

m_edu2_survey22 <- rbind(mA_edu2_survey22,mB_edu2_survey22)
m_edu2_survey22$dv <- factor(m_edu2_survey22$dv,levels=unique(m_edu2_survey22$dv))
m_edu2_survey22$lb <- factor(c("I & II","III","IV","V"), levels=c("I & II","III","IV","V"))
m_edu2_survey22

#+ fig.width=7, fig.height=4
## Figure
ggplot(m_edu2_survey22) + 
  geom_vline(aes(xintercept=1975), linetype=2, color="gray80") + 
  geom_vline(aes(xintercept=1990), linetype=1, color="gray80") + 
  geom_vline(aes(xintercept=2000), linetype=1, color="gray80") + 
  geom_vline(aes(xintercept=2010), linetype=1, color="gray80") + 
  geom_hline(aes(yintercept=0), linetype=2, color="gray50") + 
  geom_errorbar(aes(x=x, ymin=lci,ymax=uci), alpha=0.8, width=2) + 
  geom_point(aes(x=x, y=est), size=1.5, color="gray10", alpha=0.8) + 
  geom_text(aes(x=x, y=-0.7, label=lb), family="serif") + 
  facet_grid(.~dv) + 
  scale_fill_manual(values=c("black","gray70")) + 
  scale_x_continuous(breaks=c(1975,1990,2000,2010),limits = c(1963,2022)) +
  # coord_cartesian(ylim=c(-0.02,0.3)) + 
  labs(x = "Year of Turning 19 (Typical Year of Entering University)",
       y = "Conditional Coefficient with 95% Confidence Interval") +
  theme_classic(base_family="serif") + 
  theme(plot.title = element_text(hjust=0.5),
        panel.background = element_rect(color="black", size=1),
        strip.text = element_text(size=11, face="bold"),
        legend.position="none")
ggsave("mx_cplot_survey22_ologit.pdf", width=7, height=4)


#'
#' # Result Table (Table A6)
#'

screenreg(list(mA_icpp09,mB_icpp09,mA_survey22,mB_survey22), 
          digits=3, include.ci=FALSE, single.row=TRUE,
          custom.header = list("ICPP 2009"=1:2, "2022 Survey"=3:4),
          custom.model.names = c("Admit More Foreigners","Integrate Foreign Residents",
                                 "Admit More Foreigners","Integrate Foreign Residents"),
          custom.coef.map = cfmap, 
          custom.note = "%stars. Robust standard errors in parentheses. \nICPP 2009 models use standard errors clustered by residential area.")

texreg(list(mA_icpp09,mB_icpp09,mA_survey22,mB_survey22), 
       digits=3, include.ci=FALSE, single.row=TRUE,
       custom.header = list("ICPP 2009"=1:2, "2022 Survey"=3:4),
       custom.model.names = c("Admit More Foreigners","Integrate Foreign Residents",
                              "Admit More Foreigners","Integrate Foreign Residents"),
       custom.coef.map = cfmap, 
       custom.note = "%stars. Robust standard errors in parentheses (clustered by residential area in ICPP2009).",
       booktabs = TRUE, dcolumn = TRUE, use.packages = FALSE, table=FALSE, fontsize = "tiny",
       file = "mtab_ologit.tex")
